Rotating three-dimensional solitons in Bose Einstein condensates with gravity-like 

attractive nonlocal interaction 
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We study formation of rotating three-dimensional high-order solitons (azimuthons) in Bose Ein- 
stein condensate with attractive nonlocal nonlinear interaction. In particular, we demonstrate for- 
mation of toroidal rotating solitons and investigate their stability. We show that variational methods 
allow a very good approximation of such solutions and predict accurately the soliton rotation fre- 
quency. We also find that these rotating localized structures are very robust and persist even if 
the initial condensate conditions are rather far from the exact soliton solutions. Furthermore, the 
presence of repulsive contact interaction does not prevent the existence of those solutions, but allows 
to control their rotation. We conjecture that self-trapped azimuthons are generic for condensates 
with attractive nonlocal interaction. 

PACS numbers: 42.65.Tg, 42.65. Sf, 42.70.Df, 03.75.Lm 



I. INTRODUCTION 

Studies of Bose Einstein condensates (BEC) belongs to 
one of the fastest developing research directions. The ma- 
jor theoretical progress in this area has been stimulated 
by the fast experimental advances which enables to inves- 
tigate subtle phenomena of fundamental nature [l|, Q. 
In the semiclassical approach the spatial and tempo- 
ral evolution of the condensates wave function is com- 
monly described by the Gross Pitaevskii equation [3] 
which reflects the interplay between kinetic energy of 
the condensate and the nonlinearity originating from 
the interaction potential leading, among others, to the 
formation of localized structures, bright and dark soli- 
tons 0, [H]. So far the main theoretical and experi- 
mental efforts have been concentrating on condensates 
with contact (or hard-sphere) bosonic interaction which, 
in case of attraction, may lead to collapse-like dynam- 
ics. Recently, also systems exhibiting a nonlocal, long- 
range dipolar interaction jf| have attracted a significant 
attention. This interest has been stimulated by suc- 
cessful condensation of Chromium atoms which exhibit 
an appreciable magnetic dipole moment @, H @. The 
presence of spatially nonlocal nonlinear interaction and, 
at the same time, the ability to control externally the 
character of local (contact) interactions via the Fesh- 
bach resonance techniques offer the unique opportunity 
to study the effect of nonlocality on the dynamics, sta- 
bility and interaction of bright and dark matter wave 
solitons ficl [Til [T^ . [T3| . The enhanced stability of local- 
ized structures including fundamental, vortex and rotat- 
ing solitons in nonlocal nonlinear media (not necessarily, 
BEC) has been alreadypointed out in a number of the- 
oretical works 0, El El EE [3 El, H E3- In partic- 
ular, stable toroidal solitons were presented in (igL [20| . 
However, since the dipole-dipole interaction is spatially 
anisotropic, an additional trapping potential or a com- 



bination of attractive two-particle and repulsive three- 
particle interaction were necessary. Various trapping 
arrangements have been proposed to minimize or com- 
pletely eliminate this anisotropy. In particular, O'Dell 
at al. [22] have recently suggested to use a series of tri- 
ads of orthogonally polarized laser beams illuminating 
cloud of cold atoms along three orthogonal axes so that 
the angular dependence of the dipole-dipole nonlinear 
term is averaged out. The resulting nonlocal interac- 
tion potential becomes effectively isotropic of the form 
1/r. It has been already shown by Turitsyn [13] that 
a purely attractive "gravitational" (or Coulomb) inter- 
action potential prevents collapse of nonlinear localized 
waves and gives rise to the formation of localized states - 
bright solitons which could be supported without neces- 
sity of using the external trapping potential. If realized 
experimentally such trapping geometry would enable to 
study effects akin to gravitational interaction. Few re- 
cent works have been dealing with this "gravitational" 
model of condensate looking, among others, at the sta- 
bility of localized structures such as fundamental solitons 
and two-dimensional vortices [2H.l25l.l26l.l27j. 

In this paper we study formation of three-dimensional 
high-order solitons in BEC with gravity-like attractive 
nonlocal nonlinear potential. In particular, we demon- 
strate formation of vortex tororidal solitons (solitons) 
and investigate their stability. We show that such BEC 
supports robust localized structures even if the initial 
conditions are rather far from the exact soliton solutions. 
Furthermore, we also demonstrate that the presence of 
repulsive contact interaction does not prevent the exis- 
tence of those solutions, but allows to control their rota- 
tion. 

The paper is organized as follows. In Sec. [ill we intro- 
duce briefly a scaled nonlocal Gross-Pitaevskii equation 
(GPE). We discuss two different response functions, the 
above long range 1/r response and the so-called Gaus- 
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sian response yielding a much shorter interaction range. 
In Sec.HHwe recall general properties of rotating soliton 
solutions (azimuthons), which are then approximated in 
Sec. HVI bv means of a variational approach. Those vari- 
ational approximations allow us to predict the rotation 
frequency of the azimuthons which are then confronted 
with results from rigorous numerical simulations. Finally, 
self-trapped higher order three-dimensional rotating soli- 
tons are presented in Sec. [Vj and we show that such a 
nonlocal BEC's support robust localized structures. 



II. MODEL 

We consider a Bose-Einstein atomic condensate with 
the isotropic interatomic potential consisting of both, re- 
pulsive contact as well as attractive long-range nonlocal 
interaction contributions. Following O'Dell et.al [12], an 
attractive long-range interaction of "gravitational" form 
can be induced by triads of frequency detuned laser 
beams resulting in the following dimensionless Gross- 
Pitaevskii equation (GPE) for the condensate wave func- 
tion xjj (r, t): 



dtip = iAip + iOip 



e(r,t) 



(la) 
(lb) 



The nonlinear response 9 consists of both local and non- 
local contribution. Interestingly, for the "gravitational" 
nonlocal interaction contains no additional parame- 
ter (see also Appendix [Aj . The ratio between local and 
nonlocal term is solely determined by the form of the 
wavefunction ip. We will see later (Sec. IIV|1 that for very 
broad solitons the local contact interaction ~ 2 be- 
comes negligible. 

In this paper we will also consider a second, different 
nonlocal model, the so-called Gaussian model of nonlo- 
cality. Despite the fact that it is not motivated by a cer- 
tain physical system, it serves as a popular toy model 
for the general class of nonlocal Schrodinger (Gross- 
Pitaevskii) equations in one and two dimensional prob- 
lems [3,[2l,[23,[33,[3H. Here, we will extend this classi- 
cal model to three transverse dimensions, and moreover 
allow an additional local repulsive term similar to the 
previous case, and introduce 



9(r,i) 



hMr',t)|V 



d A r' -5\^ 



(2) 

The additional parameter S is necessary here to keep 
track of one of the two degrees of freedom of the Gaus- 
sian response, i.e. amplitude or width, which cannot be 
scaled out (see Appendix [Bj . The value of S determines 
the relative strength of the local repulsive term. Note 
that compared to the above "gravitational response, the 
interaction range of the Gaussian nonlocal response is 
significantly shorter due to its rapid decay for r — > oo. 



As far as stability of localized states is concerned, Tu- 
ritsyn [23] showed that the ground state of the nonlocal 
Schrodinger equation with a purely attractive 1/r ker- 
nel is stable (collapse arrest) using Lyapunoff's method. 
A rather general estimate for non-negative responsefunc- 
tions has been found in [H] for arbitrary dimensions. 
Bang at al. [H[ showed, using the same method, that 
for systems with arbitrarily shaped, nonsingular response 
functions with positive definite Fourier spectrum, col- 
lapse cannot occur. Obviously the stability of the ground 
state is only a necessary but not sufficient condition for 
the stability of rotating higher-order states, which we will 
investigate in the following by means of numerical simula- 
tions. In [23], linear and global (modulational) stability 
under small perturbations of solutions of the Hartree- 
equation was shown. 



III. ROTATING SOLITONS 

It has been shown earlier that azimuthons, i.e. multi- 
peak solitons with angular phase ramp exhibit con- 
stant angular rotation and hence can be represented by 
straightforward generalization of the usual (nonrotating) 
soliton ansatz by including an additional parameter, the 
angular frequency f2 [3~il . [3a ] . We write 



ip(r, z, 4>,t) — U (r, z, <f> — fit) 



JEt 



(3) 



where U is the complex amplitude and E is the normal- 
ized chemical potential, r — \J x 2 + y 2 and 4> denotes the 
azimuthal angle in the plane (x, y). It can be shown, that 
by inserting the above function into the nonlocal GPE 
(U) one can derive the formal relation for the rotation 
frequency [13, El 



n = - 



IL - I'M + XL — X'M 
L 2 - MM' 1 



(4) 



where the functional M, M', X, X', L, I, I' represent the 
following integrals over the stationary amplitude profiles 
of the azimuthons 



M 
L 
I 
X 

M' 
I' 

X' 



\U\ 2 d 3 r, 
U*d v Ud 3 r, 



U*AUd 3 r 



J 6(r)|C/(r)| 2 d 3 r, 
\d v U\ 2 d 3 r, 
d v U*AUd 3 r, 
i l e(r)U{d v U*)d 3 r. 



(5a) 
(5b) 
(5c) 
(5d) 
(5e) 
(5f) 
(5g) 
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The first two conserved functionals (M) and (L) have 
straightforward physical meanings of "mass" or "number 
of particles" and "angular momentum". In the next Sec- 
tion, we will compute approximate azimuthon solutions 
and their rotation frequency employing a certain ansatz 
for the stationary amplitude profile U. 



IV. VARIATIONAL APPROACH 



and the energy E is given by 
15 



E = 



2tt (49p 4 + 86p 2 + 49) - 15p ~$ ~ 15 



4ct 2 tt (49p 4 + 86p 2 + 49ct 2 tt) - 135 - 90p 2 - 135p 4 

(10) 

Because of the difference in the denominator, the local- 
ized solution (with finite amplitude) exists only if its 
width is greater than the critical value a cr (p) , 



In order to get some insight into possible localized 
states of the Gross-Pitaevskii equation we resort first to 
the so called Lagrangian (or variational) approach [ST! - 
It is easy to show that Eq. ([J) can be derived from the 
following Lagrangian density: 

C := % - («f - i>*d t iP) + | V^| 2 - ~ 6 (r, t) . (6) 

It has been shown before that rotating solitons or 'az- 
imuthons' are associated with nontrivial phase and am- 
plitude structure [3l|. In two-dimensional optical prob- 
lems the simplest case represents the state falling be- 
tween optical vortex (ring-like pattern with 2ir angular 
phase shift) and optical dipole in the form of two out- 
of-phase intensity peaks (29l. [3l|. In three dimensions, a 
reasonable ansatz for corresponding localized solutions is 



ip (r, z, (p, t) :— At exp 



„2 I Jl 



i El 



2a 2 J ' (7) 
x [cos (ip — fit) + ip sin (<p — fit)] , 

where parameter p varies between zero and unity. For 
p = Eq. (7]) describes a dipole structure consisting of 
two out-of-phase lobes, while for p = 1 it is a three- 
dimensional vortex, i.e. toroid-like structure with zero 
in the center and azimuthal (in the (x,y) plane) phase 
ramp of 2ir. Using the ansatz Eq. (7]), one can easily find 
that 



IL - I'M = 0, 



(8) 



which shows that only the non-linear terms contribute to 
the frequency d (vide formula Eq. Q). 

After inserting the solution Eq. |(7j) into the Lagrangian 
density C, and integrating over the whole 3D space we 
obtain the Lagrangian L which is the function of varia- 
tional parameters a and A only. Looking for the extrema 
of L leads to a set of algebraic relations among the vari- 
ational variables. 



A. The "gravitational" response 

In this case, the amplitude A can be expressed as a 
function of p and cr as follows (see also Appendix ICj) 



A 



5V2 (1+P 2 



49p 4 +86p 2 +49 g 
120 na 



" 4 +6p 2 +9 4 ' 



(9) 



3p 4 + 2p 2 + 3 



3 5 

2V^V 49p 4 + 86p 2 



49 



(11) 



This threshold is an obvious consequence of competi- 
tion between nonlocal and local interaction potentials, 
because the second term in the denominator of Eq. ([9]) 
is due to the local contact interaction. While the former 
being attractive, leads to spatial localization, the latter, 
which is repulsive, tends to counteract it. For small a 
the kinetic energy term is large and can be compensated 
only if the particle density is high enough. In this regime 
the local repulsive interaction prevails over the attrac- 
tion leading to the expansion of the condensate until the 
condition for its localization (i.e. a > a cr ) is satisfied. 
The rotation frequency f2 is then given by the following 
relation 



n = A 



, a 2 pV2 (4ct 2 7t-5) 



80 



(12) 



32 



Interestingly, this expression is not sign definite, which 
means that we can expect both positive and negative 
rotation frequencies. In partic ular, t he azimuthon with 
the "stationary" width <j s = w 0.63 has no an- 

gular velocity. Again, this effect is due to competi- 
tion between nonlocal and local contribution to Cl for 
small a. The nonlocal attractive interaction leads to a 
positive contribution to O, the repulsive local interac- 
tion to a negative one. The expression for fl without 
repulsion can be obtained by the outlined variational 
procedure or by asymptotic expansion (a — > oo) up to 

O (l/cr 2 ), fl = (49^4+8^+49)^^ ■ As expected, this quan- 
tity is strictly positive. Both curves fl versus a (p = 0.7) 
with and without contact interaction are shown in the 
right panel in Fig. [TJ We can see that the repulsive local 
interaction kicks in for a < 1.5. 

As we observe in Fig.Q] the mass behaves like VE close 
to E = , since generally, M ~ A 2 a 5 , and for a — > oo, 
one finds E ~ 1/c 2 M ~ l/cr, whereas for a — > 0, one 
finds E r~j l/cr 2 , M ~ cr. The fact that the mass can 
become zero for E — > is a well-known property for 
very long range kernels, such as the Coulomb potential 
in three dimensions [33j- For shorter ranged responses 
(e.g., Gaussian response, see Fig. [2]), the mass attains its 
minimum at a finite value of E. In the limit of solely 
attractive local interaction (E — > 0, a — > oo), the mass is 
a monotonically decreasing function in E. 
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FIG. 1: (color online) The left panels show the dependency 
of the mass M (top) and the width a (bottom) on the chem- 
ical potential E. Black curves show results from the varia- 
tional approach including local repulsion, dashed blue curves 
are without contact interaction. The right panel shows the 
angular frequency as a function of a. Black dots denote 
results obtained from numerical simulations of the GPE (JTJ) - 
All plots are for p = 0.7. 



B. The Gaussian response 

Repeating each step of the previous calculations for 
the Gaussian nonlocal response, one ends up again with 
expressions for amplitude A and rotation frequency f2, 
given by 



A 



9/2 



(9p 4 +9+6p 2 )tr 13 , (4p 2 + l+p 4 )o- 11 (1+p 2 ) 2 
160 ~r on ~r s 



SF Ie 



with F rcp = {a 2 + l) 9 / 2 (9p 4 + 9 + 6p 2 )a 4 /160 and 



n 



A2 p_(a 7 -S(a 2 + iy/ 2 ) a 2 V2 



16(ct 2 + l) 7 /2 



(13) 



(14) 



As already pointed out in Sec.HU the additional parame- 
ter S is necessary due to an additional degree of freedom 
of the Gaussian response, and fixes the ratio between re- 
pulsion and attraction (see Appendix [Bj. Obviously, for 
6 = 0, the repulsive local contact interaction vanishes. 

Here, a s = ^ 5 2 ' 7 / (l - S 2 / 7 ) « 0.60 for S = 0.01. 



We observe that E ~ l/a 2 ,A~ 1/cr 2 , M ~ a ~ 1/VE 
for both small and large a. Compared to the "gravi- 
tational" response, the range of this potential is much 
shorter. Hence, when considering large a the Gaussian 
response acts more and more like a local attractive re- 
sponse and higher order solitons become unstable (see 
end of Sec. EJ • 
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FIG. 2: (color online) Same as Fig. [TJ but for the Gaussian 
response given m Eq. ©. Black curves are for S = 0.01, 
dashed blue curves without repulsion (S = 0.01). All plots 
are for p = 0.7. 



V. NUMERICAL RESULTS 

In this section, the predictions of the variational ap- 
proach will be confronted with direct numerical simula- 
tions. The approximate solitons resulting from the vari- 
ational approach will be used as an initial conditions to 
our three-dimensional code to compute their time evolu- 
tion. In general, we find stable evolution, in particular 
the characteristic shape of the initial conditions is pre- 
served. For rotating azimuthons, the angular velocities 
will be measured and compared to the ones obtained in 
the previous section. 

In Fig. [3] we illustrate the temporal evolution of three- 
dimensional solitons for the "gravitational" response, i.e., 
solutions to Eq. (I]) . This first two rows present the clas- 
sical stationary soliton solutions torus and dipole, respec- 
tively. Due to imperfections of the initial conditions ob- 
tained from the variational approach we observe slight 
oscillations upon evolution, in particular for the dipole 
solutions (second row). Those oscillations are not present 
if we use numerically exact solutions (obtained from an 
iterative solver [38]) as initial conditions (not shown). In 
the last row of Fig. [3] we show the evolution of an az- 
imuthon (p = 0.7, rj = 1). The rotation of the amplitude 
profile is clearly visible. Again we observe radial oscil- 
lations due to the imperfect initial condition, but the 
solution is robust. 

Figure [4] shows the dependency of the azimuthon ro- 
tation frequency as a function of the modulation param- 
eter p. Solid lines represent predictions from the vari- 
ational model, black dots represent rotation frequency 
obtained from numerical simulations. As expected from 
two-dimensional nonlocal models [1q|,[39J], the modulus of 
f2 increases with p. Our variational calculations predict 
that for small width a, when repulsive interaction comes 
into play, the sense of azimuthon rotation changes. In 
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FIG. 3: (color online) Dynamics of the three-dimensional sta- 
ble solitons in gravity-like BEC. Iso-surfaces of the normal- 
ized density |i/>| 2 are depicted for different evolution times, 
the interior density distribution is represented in grey-scales. 
The initial variational parameters used are a = 1 and p = 1 
(torus, iso-density surface at ji/^ 2 = 0.76) for the upper row, 
p = (dipole, iso-density surface at |i/>| 2 = 1.41) for the mid- 
dle one and finally p = 0.7 (azimuthon, iso-density surface at 
l^] 2 = 0.86). The sense of the rotation (f2 = 0.64) is indicated 
by the arrows. 



particular, we found a "stationary" width a s where the 
rotation frequency fi vanishes. Indeed, full model simu- 
lations confirm this property, since the first row in Fig. [5] 
and HI a) show a very slow rotation with opposite orien- 
tation, so that the numerical stationary width is between 
0.6 . . . 0.61. Hence, we propose that tuning the strength 
of contact interaction in experiments allows to control 
the azimuthon rotation. 

Furthermore, we observe that very narrow azimuthons 
(a — > a cr ) have negative f2 and rotate very fast (see 
Fig. [I]). This may be interesting for potential exper- 
iments, since the duration of BEC experiments is re- 
stricted to typically several hundreds of milliseconds. 
However, for azimuthons very close to cr cr the ansatz 
function [7] becomes less appropriate and using variational 
initial conditions leads to very strong oscillations upon 
evolution, up to the point where is is no longer possible 
to identify properly the rotation frequency f2. 

Concerning the Gaussian nonlocal response, we find 
very similar evolution scenarios. Results shown in the left 
panel in Fig. [6] for the Gaussian response underline the 
observations from above, in particular we also find non- 
rotating azimuthons at a = a s . However, there are some 
important differences. First, it seems that our ansatz 
Eq. [7| is better suited for the Gaussian response, the ra- 
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FIG. 4: (color online) Azimuthon rotation frequency f2 vs 
modulation parameter p in gravity-like BEC, for a = 0.6 ss a 3 
(left panel) and a = 3 (right panel). Black curves show re- 
sults from the variational approach including local repulsion, 
dashed blue curves are without contact interaction. Black 
dots denote results obtained from numerical simulations of 
the GPE p. 






FIG. 5: (color online) The upper row shows iso-density sur- 
faces at |^j 2 = 7.63 for the very slow rotating (f2 m 0) az- 
imuthon with p = 0.7 and a = 0.61. The lower row shows a 
fast counter-rotating (fl = —2.24) azimuthon with p — 0.7, 
a = 0.5 and iso-density surface at \ip\ 2 = 32. Same plot style 
as in Fig. [3] 



dial oscillations we observed with the "gravitational" are 
still present, but much weaker. The second difference is 
due to the fact that the Gaussian response has a much 
shorter range than the "gravitational" one. For large a 
the Gaussian kernel acts like a attractive local response. 
As a consequence, higher order solitons become unstable 
in the sense that the two humps spiral out. We observe 
unstable evolution in numerical simulations for a > 0.9 at 
p = 0.7. The right panel in Fig. [6] visualizes the cause of 
this instability: For increasing sigma the resulting con- 
volution term 8 [Eq. Q], which is responsible for the 
self-trapping, becomes smaller in amplitude and asym- 
metric in the rotation plane, which eventually leads to 
destabilization of the azimuthon. 
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FIG. 6: (color online) Azimuthon rotation frequency vs 
modulation parameter p in BEC with solely attractive Gaus- 
sian nonlocal response for a = 0.6 » a s (left panel). Black 
curves show results from the variational approach including 
local repulsion (S = 0.01), dashed blue curves are without con- 
tact interaction (5 = 0). Black dots denote results obtained 
from numerical simulations (Eq. [2j 5 = 0.01). The right panel 
shows profiles of the convolution term O [Eq. ([2])] for a — 0.6 
and a = 1. Solid lines correspond to profiles along the major 
axis of the resulting ellipsiod (z — 0, ip = 0), dotted lines to 
those along the minor axis (z = 0, tp = 7r/2). 



VI. CONCLUSION 

We studied formation of rotating localized structures 
in Bose Einstein condensate with different nonlocal inter- 
action potentials. We successfully used variational tech- 
niques to investigate their dynamics and showed numer- 
ically that such localized structures are indeed robust 
objects which persist over long evolution times even if 
the initial conditions significantly differ from the exact 
soliton solutions. 

For rotating solitons (azimuthons) , we derived analyti- 
cal expressions for the angular velocity, in excellent agree- 
ment with rigorous three-dimensional numerical simula- 
tions. Furthermore, we show that it is possible to control 
the rotation frequency by tuning the local contact inter- 
action, which is routinely possible by Feshbach resonance 
techniques. In particular, we can change the sense of ro- 
tation, and we can find non-rotating azimuthons. We also 
identify parameter regions with particularly fast rotation, 
which may be important for potential experimental ob- 
servation of such solutions. 

By using different nonlocal kernel functions we showed 
that rotating soliton solutions are generic structures in 
nonlocal GPE's. Hence, we conjecture that the phenom- 
ena observed in this paper are rather universal and apply 
for a general class of attractive nonlocal interaction po- 
tentials. 



APPENDIX A: NORMALIZATION OF THE 
"GRAVITATIONAL" MODEL 

We consider a Bose-Einstein atomic condensate with 
isotropic interatomic potential consisting with both, re- 
pulsive contact as well as attractive long-range nonlo- 
cal interaction contributions. Following O'Dell et.al 



the attractive long-range interaction which is electro- 
magnetically induced by the triads of frequency detuned 
laser beams with the intensity I can be presented in the 
"gravitational" form 



U Gr 



11 Iq 2 a 2 1 
47r ceo f 



(Al) 



Here, q is the modulus of the wave vector, a the isotropic, 
dynamic polarizability of the atoms, c the light velocity 
and eo the permittivity of the free space. Then the com- 
plete two-body interaction potential is given by 



V(r) = S(r) - -, 

m r 



(A2) 



where the first term comes from the contact s-wave scat- 
tering, a is the scattering length, and m is the atomic 
mass. The potential can only be written in this form, if 
the mean kinetic energy per particle dominates the u/f- 
term, so that the short-range hard-sphere scattering is 
not affected. This is fulfilled, if a B <C A B <C a, where 

,2 

a B '■= rrr: is the Bohr radius, associated with the inter- 
action, and Ab is the de Broglie wavelength. 

The temporal and spatial dynamics of the condensate 
wave function ip (f,i) is then governed by the the follow- 
ing Gross-Pitaevskii equation (GPE): 



iHdxip + - — A}tp + u 
2m 



Using the normalization 



|r - r', 

Anah 2 
m 



(A3) 



2 _ 

ip = 0. 



ib := J8iraR 2 i> =: — $ 



um _ 1 
r =: — 

4:Trah 2 R c 



h - 1 ~ 

-t =: —t, 

2mR 2 T c ' 



(A4a) 
(A4b) 
(A4c) 



one ends up with the dimensionless GPE {TJ. The ac- 
tual values for the scaling parameters are R c — 19/im, 
m = 3.8 • l(T 26 kg, a = 3nm, u = 2 • 10- 13 eVnm, 
T c = 0.25s which corresponds to typical experimental 
conditions [22]. Then one finds that for the condensate 
consisting of, say, N = 10000 atoms the above normal- 
ization gives 



M = / |^| d 3 r = 



N 



42. 



(A5) 
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APPENDIX B: NORMALIZATION FOR THE 
GAUSSIAN MODEL 



APPENDIX C: CONVOLUTION 



For the Gaussian response, one has to start from the 
equation 



ihdfip H Aib 

2m 

4>(r',t) 



4nah 2 



+a r 



(Bl) 



r d d ?V = 



Here, the degree of nonlocality (i.e. <Jr), that is fixed 
for the 1/r-response, and the amplitude Ar of the the 
response function can be chosen. By using the normal- 
ization 



r 

t 
6 



ft2 3 / 4 



1 

4 a 2 TT 3 / 4 ^ARam 
1 . 



1 - 



2a m ~ 



1 - 

1 =: —t 

ah 2 V2 



o^^pK-mAR 



(B2a) 
(B2b) 
(B2c) 
(B2d) 



one ends up with Eq. [2j that has the additional degree of 
freedom S. 



The convolution term in the Lagrangian ([6]) can be 
calculated analytically by, for instance, re-writing the 
integrand in terms of spherical harmonics Yy . Then, 
|^(r)| 2 = ^VijYij, where denote the coefficients of 
the spherical harmonics, and the convolution integral can 
be easily calculated leading to the following result: 



d 3 r'd 3 r 



o Jo 



AirrA r exp 



x U 2 + ^lo(7) 2 + ^i ±2 (j^ I r r \\r'r-dr 



OO POO 



Jr 



Anr 1 A r 2 exp 



x(yoo + ^lo(^) 2 + ^ 2 2 ±2 (^) 2 ) 



r' dr'r dr 
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= 7777; <7 7T 2 A — =. 

240 y/2 
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